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ABSTRACT 


A finite-difference  solution  technique  has  been  developed 
for  subsonic  two-dimensional  inviscid  flow  past  lifting  airfoils. 
This  work  is  an  adaptation  of  the  method  used  by  Sells  (1967). 
The  full  governing  equations  of  compressible  flow  are  written  in 
terms  of  a translated  velocity  potential  which  is  continuous 
throughout  the  flow  field.  This  simplifies  solutions  for  bluff 
airfoils  (no  Kutta  condition)  where  both  angle  of  attack  and 
lift  coefficient  are  specified.  The  computational  plane  is  the 
interior  of  a unit  circle  obtained  by  mapping  the  flow  field  into 
the  interior  of  the  circle.  A line  overrelaxation  matrix  method 
is  used  for  solution  of  the  partial  differential  equation  which  in 
the  iteration  scheme  is  coupled  with  an  algebraic  equation.  The 
numerical  procedure  is  accurate  and  well  behaved  for  all  sub- 
sonic flow  conditions. 


ADMINISTRATIVE  INFORMATION 

The  work  presented  herein  was  conducted  for  the  Office  of  Naval  Research  (Code  460) 
as  Project  Order  2-0127,  NR  212-204  and  was  accomplished  in  the  time  period  July  1971  to 
April  1973. 

The  material  was  previously  submitted  to  the  University  of  Maryland  in  partial  fulfillment 
of  requirements  for  the  degree  of  Master  of  Science,  Aerospace  Engineering.  Thus,  in  some 
details  the  report  deviates  from  the  traditional  format  of  the  Naval  Ship  Research  and 
Development  Center  (now  David  W.  Taylor  Naval  Ship  Research  and  Development  Center). 
Preparation  of  this  report  was  funded  under  Work  Unit  1-1619-111. 

INTRODUCTION 

An  analysis  of  the  inviscid  flow  of  airfoil  sections  is  of  considerable  value  in  designing 
new  sections  and  in  understanding  the  performance  of  existing  airfoils.  Although  viscous 
effects  influence  airfoil  characteristics,  the  distribution  of  inviscid  pressure  generally  agrees 
closely  with  experimental  results,  especially  on  the  forward  60  to  80  percent  of  the  airfoil. 
Boundary  layer  control  techniques,  which  are  becoming  increasingly  common,  reduce  the 
effects  of  boundary  layers  and  thus  make  the  inviscid  condition  even  more  applicable.  Despite 
the  usefulness  of  such  solutions,  no  method  has  been  available  for  an  exact  solution  of  the 
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governing  equations  of  subsonic  compressible  flow.*  The  difficulty  arises  from  the  nonlinear 
nature  of  the  equations. 

Incompressible  solutions  for  arbitrary  airfoils  have  been  available  for  many  years.  Until 
very  recently,  however,  compressible  flow  solution  methods  have  been  limited  to  the  application 
of  so-called  compressibility  correction  factors  to  the  incompressible  solution.  The  best  known 
of  these  are  the  Prandtl-Glauert  and  Karman-Tsien  factors  that  are  derived  from  linearization 
of  the  full  equations.  It  is  generally  recognized  that  they  are  adequate  for  flow  conditions 
well  below  those  that  produce  sonic  velocity.  Unfortunately,  accuracy  is  lost  as  local  sonic 
velocity  is  approached  - just  where  compressibility  effects  are  highest.  The  error  is  greatest 
near  the  leading  edge  where  the  perturbation  velocities  are  largest  and  at  any  point  where 
local  sonic  velocity  is  approached.  This  deficiency  is  severe  enough  to  prohibit  even  an 
approximate  determination  of  the  critical  Mach  number  of  airfoils  with  rapid  flow  acceleration 
around  the  leading  edge  as  on  blunt-nosed,  “peaky”  type  sections. 

In  1967  Sells  (Reference  1)  presented  a practical  numerical  scheme  for  solving  the  sub- 
sonic flow  equations.  The  stream  function  was  used  as  the  unknown  variable.  Unfortunately, 
the  behavior  of  the  stream  function  at  sonic  velocity  is  such  that  increasingly  severe  numerical 
problems  are  encountered  as  sonic  velocity  is  approached.  Solutions  are  attainable  up  to  a 
maximum  local  Mach  number  of  about  0.95  and  0.98.  Even  for  this,  computational  time 
becomes  excessive  and  error  develops  because  of  the  polynomial  curve  fit  required  when  this 
approach  is  used. 

This  study  originated  as  an  attempt  to  avoid  the  limitations  imposed  by  the  stream 
function  while  retaining  some  of  the  desirable  concepts  used  by  Sells.  It  has  led  to  a practical, 
well-behaved  numerical  technique  for  the  exact  solution  of  two-dimensional  subsonic  flow. 

This  success  resulted  from  using  the  velocity  potential  rather  than  the  stream  function. 


*Shortly  after  this  project  was  completed,  it  was  learned  that  Garabedian  and  Korn  of  New  York  University 
were  publishing  a numerical  technique  for  the  compressible  flow  equations.  Time-dependent  solutions  were 
demonstrated  several  years  ago  by  Yoshihara  and  Magnus  and  others. 


APPROACH 


BASIC  EQUATIONS* 


Flow  conditions  will  be  considered  as  steady,  inviscid,  isentropic,  and  irrotational.  The 
continuity  equation  can  be  expressed  as 

V • (p  • u)  = 0 (1) 

with  the  irrotationality  condition  given  by 

V x u = 0 (2) 

The  Bernoulli  integral  of  the  energy  equation  is 

a2  1 _ 

7 + - u2  = constant 

7-1  2 


With 


the  isentropic  relationship  of 


Poo  = 1 , I Uqo  I = 1 , 3oo  = 1 /M 


ii  = _L 

2 Poo  P \Pooj 


can  be  used  to  give 


a2  = pT-i/M2 


The  Bernoulli  equation  then  becomes 


PJ~l  ^ * —,2 I ,1 


+ ^ iu  r = 


(7- l)M2  2 (7  - 1 )M2  2 


(3) 


ADVANTAGES  OF  THE  VELOCITY  POTENTIAL 

The  above  equations  imply  the  existence  of  both  a velocity  potential  and  a stream 
function.  Equations  could  be  written  with  either  of  these  two  as  the  unknown  variable. 


*The  notation  of  Reference  1 is  generally  retained  here  for  convenient  comparison  with  the  corresponding 
development  for  the  stream  function. 
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Figure  1 — Density  as  a Function 
of  Mass  Flow 


Figure  2 - Density  as  a Function 
of  Velocity 


As  mentioned  previously,  the  stream  function  is  ill-behaved  as  local  sonic  velocity  is  approached. 
The  nature  of  this  problem  can  be  seen  in  Figure  1.  The  mass  flow  parameter  (which  is 
associated  with  stream  tubes)  versus  density  relationship  is  multivalued  with  a saddle  point  at 
the  speed  of  sound.  Beginning  at  about  M = 0.85,  numerical  instabilities  are  encountered 
because  of  the  increasingly  large  changes  in  density  with  small  changes  in  mass  flow.  Figure  2 
illustrates  the  behavior  of  the  velocity  potential.  This  relationship  is  single  valued  and  well- 
behaved.  The  potential  is  therefore  to  be  preferred;  however,  it  has  a characteristic  which 
must  be  dealt  with  very  carefully  in  a numerical  scheme. 

There  is  a discontinuous  change  in  value  of  the  potential  across  a line  connecting  the 
rear  stagnation  point  to  infinity.  The  jump  across  this  Kutta  slice  is  proportional  to  the 
circulation.  This  can  be  readily  handled  numerically  if  the  location  of  the  rear  stagnation 
point  is  known.  From  the  Kutta  condition,  this  is  the  trailing  edge  point  on  conventional 
airfoils.  In  recent  years  there  has  been  increasing  interest  in  airfoils  with  bluff  trailing  edges 
(no  Kutta  condition).  For  these  sections,  circulation  is  established  and  controlled  by  boundary 
layer  control  techniques  such  as  tangential  slot  blowing  over  the  rounded  trailing  edge. 

Analysis  of  the  inviscid  flow  of  such  airfoils  requires  the  specification  of  both  angle  of  attack 
and  lift  coefficient.  The  rear  stagnation  point  is  determined  as  part  of  the  solution;  hence  the 
location  of  the  Kutta  slice  is  not  known  a priori.  This  difficulty  will  be  handled  quite  simply 
by  introducing  a translated  potential  which  is  continuous  throughout  the  flow  field. 

NUMERICAL  SOLUTION  TECHNIQUE 

To  minimize  the  complexity  of  the  finite  difference  solution  that  must  be  performed, 
the  flow  equations  are  written  in  the  form  of  a partial  differential  equation  (PDE)  coupled 


with  an  algebraic  expression.  In  the  iteration  scheme,  line  overrelaxation  matrix  methods  are 
used  to  solve  the  continuity  PDE  for  0 while  the  density  field  is  held  constant.  Then  the 
density  at  all  grid  points  is  recalculated  from  the  Bernoulli  equation  by  using  the  new  value 
of  0.  This  process  is  repeated  until  successive  iterations  produce  insignificant  changes  in  the 
unknown  variables. 

COORDINATE  SYSTEM 

The  transformed  coordinate  system  as  used  by  Sells  will  also  be  applied  here.  This 
system  employs  conformal  mapping  methods  such  as  that  of  Theodorsen  to  map  the  airfoil 
into  a unit  circle.  The  exterior  (flow)  field  is  then  transformed  by  a 1/R  function  into  the 
interior  of  the  circle  and  this  has  several  advantages.  Infinity  now  appears  at  the  circle  center 
(r  = 0),  and  that  enables  the  entire  infinite  flow  field  to  be  included  in  the  computations.  A 
constant  Ar,  A 6 grid  mesh  in  the  working  plane  (Figure  3)  provides  convenient  representation 
of  the  finite  difference  expressions  with  the  airfoil  surface  always  forming  one  side  of  an 
essentially  rectangular  computational  cell.  The  most  important  feature  is  that  this  simple  grid 
results  in  an  inherent  mesh  refinement  in  the  physical  plane.  This  is  readily  apparent  from 
Figure  3.  The  grid  mesh  is  of  higher  density  near  the  surface  and  the  leading  and  trailing  edges 
just  where  the  flow  gradients  are  largest.  This  ensures  that  the  nose  region  is  calculated 
properly.  This  has  been  an  area  of  deficiency  in  many  other  attempts  at  fluid  flow 
computations. 


GRID  IN  THE  COMPUTATIONAL  PLANE  GRID  IN  THE  PHYSICAL  PLANE 


(From  Sells,  Ref.  1) 
Figure  3 — Grid  System 
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DEVELOPMENT  OF  EQUATIONS 


EQUATIONS  IN  TRANSFORMED  PLANE 


The  governing  equations  are  now  to  be  transformed  into  the  unit  circle  coordinates  r,  9 
by  using  orthogonal  curvilinear  coordinate  system  principles.  Curvilinear  scale  factors  are 
represented  by  hj  and  h2.  The  element  of  length  ds  in  the  Z-plane  is  given  by 

ds2  = | dz  |2  = hj2  dr2  + h22  dd2 


I dz  |2  = 


Idol2  = B2  (dr2  + r2  d02) 


(4) 


from  which  we  find 


h1  = B,  h2  = rB 


(5) 


where  B(r,  9)  is  the  transform  modulus  | dz/da  | between  the  physical  and  computational 
planes  as  determined  by  the  airfoil  mapping  process.  The  continuity  equation  (1)  becomes 


|;(ph2Ur)  + ^ (phi  Ud)  = 0 


(6) 


From  the  irrotationality  condition  (2) 


9 9 

a^(h2u0)-^  (hlUr)  = 0 


(7) 


there  is  a function  0(r,  0)  such  that 


00  90 

39  h2u0'  hl  ur 


(8) 


Substituting  (8)  into  (6), 


0 ( h2 


0 ( hl 


0r  \P  h.  ^ + dd  \P  h. 


= 0 


(9) 


or,  from  (5),  the  continuity  equation  is 


9_  / 90 \ , J.  (p  90\  = n 
0r  \Pf  0r/  00  Vr  00/  ° 


(10) 
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The  Bernoulli  equation 


P=|  l+YM»(l-u2) 

remains  unchanged  except  for  the  expression  for  velocity  obtained  from  (5)  and  (8): 


(ID 


I n |2  - „ 2 4.  „ 2 _ _L 


ur  + V = 


'30  2 +J_  90 
9r  j.2  9 0_ 


(12) 


BOUNDARY  CONDITIONS 

At  the  surface,  appropriate  boundary  conditions  are  no  flow  normal  to  the  surface 

ur  = 0 (13) 

and  no  flow  around  the  sharp  trailing  edge  (if  there  is  one), 

0q  = 0 for  0 = 0 (14) 

The  transformation  procedure  is  assumed  to  map  the  trailing  edge  into  the  point  0 = 0.  This 
boundary  condition  is  not  applied  for  bluff  airfoils. 

Conditions  at  infinity  correspond  to  free-stream  values  with 

P — Poo?  ^ ^oo 

The  infinity  condition  of  free-stream  velocity  cannot  be  imposed  directly;  what  is  needed 
is  an  expression  for  0.  For  a first  examination  of  0 at  infinity,  we  assume  that  at  large 
distances  | Z | in  the  physical  plane,  the  density  is  essentially  the  free-stream  value  and  there- 
fore can  be  considered  constant.  The  flow  is  then  essentially  incompressible  so  that  we  can 
examine  the  complex  potential  £2  for  steady  flow  at  incidence  a with  circulation  represented 
by  E.  For  large  | Z | (see  any  standard  text  on  incompressible  flow) 

' £2  = 0 + i0  ~ Ze“*a  - i E In  z 

With  Z = Re1^,  the  velocity  potential  can  be  expressed  as 

0 = R cos  (v?  - a)  + v?E  (15) 

and  with  transformation  into  the  circle  by  R ~ 1/r,  = -0 

0(0,  0)  ~ -p  cos  (0  +a)  - 0E  (16) 

Thus,  conditions  at  infinity  are  represented  by  a dipole  and  a vortex. 
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TRANSLATED  POTENTIAL 


The  previous  expression  for  0 at  infinity  reveals  a singularity  at  r = 0 and  a jump  in  0 at 
0 = 27 r.  The  singularity  and  the  jump  can  be  removed  by  defining  a translated  potential 

$(r,  0)  = 0 - — cos  (0  + a)  + 0E  (17) 

r 

This  function  is  finite  and  continuous  everywhere.  The  numerical  scheme  will  solve  for  d>, 
and  if  needed,  the  value  of  0 itself  can  of  course  be  found  from 

0(r,0)  = $ + - cos(0  + a)-0E  (18) 

r 

The  boundary  condition  at  infinity  for  $>  will  be  found  by  first  taking  a closer  look  at  the 
asymptotic  behavior  of  0. 


BOUNDARY  CONDITION  FOR  THE  TRANSLATED  POTENTIAL 


Expanding  0 and  p about  r = 0 to  one  more  term,  we  have 

0 = — cos  (0  +a)  - 0E  + L(0)  + 0(r) 
p = 1 + rK(0)  + 0(r2) 


(19) 

(20) 


In  addition  it  should  be  noted  that  for  large  I Z | , the  nature  of  the  transformation  employed 
is  such  that  (Reference  1)  - 


B = 


dz 

da 


\ [1  + 0(r2)] 
r2 


Substituting  the  above  in  the  continuity  equation  (10),  terms  0(l/r2)  cancel  and  terms  0(1  /r) 
give 

^ [L'(0)-K(0)sin(0  +«)]  = 0 (21) 

so  that 

L'(0)  - K(0)  sin  (0  + a)  = -jSj  (22) 

Similarly  with  the  expansion  for  0 and  p substituted  into  the  Bernoulli  equation,  terms 
0(r)  give 
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K(0) 

M2 


+ E sin  (0  +a)  - L'  sin  (0  + cO  - 0 


(23) 


Equations  (22)  and  (23)  can  be  solved  for  K and  L' 

(E  +0j)  sin  (6  +<x) 


K(0)  = 


L'(fl)  = 


sin2  (0  + a)  - ( — ) 

\M2/ 

E sin2  (0  + a)  +/3, 

\M2/ 


(24) 


sin2  (0  + a) 


(25) 


Integration  of  (25)  leads  to 


-(E+/V  _i 
L(0)  = — ; ==  tan 


y I - M2 


yj  1 - M2  tan  (0  + a )J  + E(0  + q)  + j3 


(26) 


The  additive  constant  {3  has  no  physical  significance  in  the  potential  and  can  be  dropped 
without  loss  of  generality  . From  (19),  (26),  and  ( 17).  the  boundary  value  ot  becomes 


<t>(0,  0)  = 


1 E '*'01  ) -tr 


y i - m; 


tan 


y 1 - M2  tan  (0  +<*)J  +E(0  + 


a) 


(27) 


This  asymptotic  value  is  a function  of  the  direction  from  which  r = 0 is  approached.  The 
constant  (31  will  be  determined  from  the  consideration  that  be  continuous  with  no  jump  at 
the  rear  stagnation  point  0S  (2n  for  conventional  airfoils).  Observe  that  1.(0)  -*■  0 as  M -*■  0 
is  also  required. 

CIRCULATION  CONSIDERATIONS 

From  the  definition  of  circulation  we  have 

*0C  + 2 7T 


ds  = J 


if,'™ 


-1 


0c  + 27 r 


90 

90 


d0 


(28) 


where  0S  is  the  location  of  the  rear  stagnation  point. 

Obtaining  0#  from  (18).  we  then  have 

r = - 2tt  E + <J>  (0g  + 2tt)  - 4>  (0g) 


(29) 
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and  by  using  (27),  one  obtains 


-(E+/J.) 

r = ■ 2tt 

x/l-M2 


The  value  of  remains  to  be  found.  It  is  desired  that  $ be  a continuous  function 
which  implies  that 


/!?*■/ 


0S  + 27T  g^ 


dd  = 0 


From  this  requirement 


^ = e (y  i - m2  - i) 


r = -27rEorE  = - — 

2it 

Finally,  with  these  constants  we  obtain  the  boundary  value  of  $ 

$(r->  0,  0)  = jtan  [\/  1 - M2  tan  (6  +a)J  - (6  +a)| 

The  relationship  between  circulation  and  the  lift  coefficient  Cg  is  given  by 

/ 1 2 
L — Poo  Uqq  P — Cg  C ~ Poo  Uoo 

where  c'  is  the  chord  of  the  airfoil  mapped  to  the  unit  circle.  The  result  is  that 


r _ ZZi! 

4 o' 


or  for  bluff  airfoils  where  Cg  is  a specified  condition, 

-C0c' 
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FINAL  EQUATIONS 


Equations  in  terms  of  4>  are  readily  obtainable  by  substituting  (18)  into  (10).  With 


4>r  cos  (0  + a)  + $>r 


= - — sin  (0  + a)  - E + $a 
r u 


0rr  = — cos  (0  + a)  + <£ 
r3 


rr 


= - — cos  (0  + a)  + $qq 


we  obtain 


sin  (0  + a)  + E 


£r2  4>r  - cos  (0  + o!)J  pr  - ^ sir 

' 00  +P  [f2  $rr  + r + ^0  J = 0 


(38) 


This  is  the  PDE  to  be  solved  by  matrix  methods  while  the  density  p is  held  constant.  The 
density  is  then  recalculated  from  the  Bernoulli  equation  (11)  by  using  the  newly  obtained 
values  of 


p = 


1 


1-1 

B2 


1 


<I> cos  (0  + a.) 

I 7 


2 r 


2 \-f]l/(7-l) 


— sin  (0  + a) 


-Er) 


Surface  boundary  conditions  become 

3>r  = cos  (0  + a) 

If  there  is  a Kutta  condition 

U0  = - g [sin  (0  + a)  + E - <J>0]  =0 

from  which  follows  the  means  by  which  the  circulation  is  established: 

E = $0  - sin  (a) 


(39) 


(40) 


(41) 


For  bluff  airfoils,  E is  calculated  directly  from  (37)  by  using  the  input  Cg  with  the  rear 
stagnation  point  location  found  as  part  of  the  solution.  Conditions  at  infinity  are  given  by 
(27). 
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NUMERICAL  METHODS 


The  complete  numerical  scheme  is  diagrammed  in  Figure  4.  The  discussion  in  the  next 
two  sections  refer  to  the  MATRIX  SOLN  block  (Figure  4). 


Figure  4 — Flow  Diagram  for  Numerical  Solution 

FINITE  DIFFERENCES 

Figure  5 is  a schematic  of  the  computational  grid.  For  each  point  we  write  the  finite- 
difference  equivalence  of  the  derivatives  in  the  continuity  equation  (38).  The  equation  is 
elliptic  for  subcritical  flow;  this  means  that  point  0 in  Figure  5 is  influenced  by  all  surrounding 
points.  The  finite  differencing  must  reflect  this  domain  of  dependence  by  using  centered 
differences.  Substitution  of  the  finite-difference  expressions  in  (38)  will  result  in  an  equation 
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UPPER  SURFACE 


LEADING 

EDGE 


Figure  5 — Finite  Difference  Computational  Grid 

for  $ at  the  point  0 in  terms  of  p and  $ at  adjacent  points.  (Details  of  this  development  are 
given  in  the  appendix.)  This  is  done  for  each  point  in  the  field  of  typically  1 800  points, 
resulting  in  a like  number  of  simultaneous  equations  for  the  values  of  <h. 

LINE  OVERRELAXATION 

It  is  not  feasible  to  solve  the  complete  set  of  equations  simultaneously.  Instead  an 
iteration  scheme  can  be  formed  to  enable  the  flow  field  to  be  solved  in  successive  parts. 
Referring  to  Figure  5,  we  can  use  standard  matrix  methods  to  solve  for  <f>  on  successive  lines  - 
either  lines  of  constant  6 or  constant  r.  Lines  of  constant  r,  that  is  concentric  circles,  were 
found  to  be  more  efficient  in  computer  time.  With  an  initial  condition  of  incompressible 
flow,  a Gauss-Seidel  matrix  method  is  used  on  each  line  starting  from  the  innermost  circle 
(nearly  free-stream  conditions)  and  working  to  the  one  adjacent  to  the  surface.  The  values  of 
<[>  newly  computed  for  one  line  are  used  on  the  line  that  follows,  and  so  on.  Surface  values 
of  $ are  found  from  the  enforcement  of  the  flow  tangency  condition  thereby  using  one-sided 
differences  (see  the  appendix). 

Line  overrelaxation  was  found  to  greatly  reduce  the  number  of  iterations  required.  (An 
iteration  consists  of  one  pass  through  the  complete  computational  grid  to  obtain  new  estimates 
of  $ followed  by  the  recalculation  of  density  at  all  points.)  Overrelaxation  accelerates 
convergence  by  using  the  difference  between  successive  iterations  to  estimate  a presumably 
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more  accurate  value.  This  is  done  on  a line-by-line  basis  and  is  expressed  by 

<t>  = <b  03  + (1  - to)  $ 

new  new  v ’ previous 

where  co  is  the  relaxation  factor.  Theoretically,  values  approaching  2.0  can  be  used  without 
instabilities  developing.  It  was  found  that  a factor  of  1.9  produced  the  fastest  convergence. 
The  number  of  iterations  for  convergence  is  shown  in  Figure  6 as  a function  of  the 
relaxation  factor. 


0 i i i 1 1 1 1 1 1 1 1 

1.0  1.2  1.4  1.6  1.8  2.0 

RELAXATION  FACTOR 


Figure  6 — Number  of  Required  Iterations  as  a Function 
of  Relaxation  Factor  for  a Typical  Case 

Radial  line  (constant  Q)  relaxation  was  also  tried.  Since  the  boundary  conditions  at  each 
end  of  the  line  are  known,  a simple  tridiagonal  matrix  solution  is  possible.  The  radial  line 
approach  was  found  to  be  much  less  satisfactory,  as  will  be  discussed  later. 

CONVERGENCE 

The  overall  numerical  scheme  has  been  diagrammed  in  Figure  4.  After  a pass  through 
the  field  to  calculate  new  values  of  <I>,  the  density  at  each  point  is  updated  by  using  the 
Bernoulli  equation.  Circulation  is  recalculated  before  the  density  calculation.  Then  another 
pass  is  made  for  <f>.  When  the  density  is  updated  at  each  point,  the  extent  of  change  from 
the  previous  value  is  checked.  When  this  difference,  or  residual,  is  below  a specified  level, 
convergence  is  considered  complete.  It  was  found  that  a density  change  of  0.25  x 10~4 
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(normalized  to  free-stream  density)  gave  sufficient  accuracy  for  practical  applications. 

Figure  7 shows  how  the  number  of  iterations  required  for  a given  residual  varies  for  a typical 
case.  Residuals  of  circulation  or  $ could  also  be  used  as  convergence  criteria;  however,  they 
were  not  tried  here. 


Figure  7 — Residual  Reduction  with  Increasing  Iteration 
for  a Typical  Case 

COMPUTATIONAL  ASPECTS 

It  was  found  that  time  per  iteration  on  the  CDC  6600  computer  can  be  expressed  by 

T (sec/iteration)  = 0.24  (£m)  x 10"3 

where  £ is  the  number  of  radial  lines  and  m the  number  of  concentric  circles  or  increments 
in  6 and  r,  respectively.  The  number  of  iterations  required  for  convergence  in  a given  case 
appeared  to  be  linearly  proportional  to  m;  thus  the  total  time  for  a solution  is  proportional 
to  £m2. 

The  mesh  size  required  to  obtain  a given  accuracy  depends  on  the  gradients  of  the 
dependent  variable.  The  grid  selected  as  the  best  compromise  between  computational  costs 
and  accuracy  was  £ = 120,  m = 15.  If  the  pressure  distribution  is  smooth  (“non-peaky”)  then 
a grid  of  120  x 10  is  acceptable.  For  peaky,  high-leading-edge  acceleration  profiles,  a grid  of 
160  x 15  is  recommended. 
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MACH  NUMBER 

Figure  8 - Variation  of  Number  of  Iterations 
with  Mach  Number  for  a Typical  Case 

The  total  number  of  iterations  is  also  a function  of  the  extent  of  compressibility  effects 
which  are  proportional  to  the  free-stream  Mach  number  and  the  peak  surface  velocity. 

Figure  8 illustrates  the  variation  of  number  of  iterations  with  Mach  number  for  a typical  case. 

The  above  discussion  applies  to  the  initial  conditions  corresponding  to  incompressible 
flow.  If,  for  example,  several  Mach  numbers  are  to  be  run  for  a given  configuration,  then  the 
converged  solution  can  be  used  as  the  starting  condition  for  the  next  Mach  number.  The 
number  of  iterations  required  is  thus  substantially  reduced.  This  technique  offers  one  way  of 
checking  the  degree  of  convergence.  This  is  done  by  computing  the  same  case  by  starting 
from  converged  solutions  below  and  above  the  free-stream  Mach  number  of  interest. 

As  mentioned  previously,  radial  line  overrelaxation  was  also  investigated.  At  first  this 
was  tried  by  starting  at  6 = 0 and  moving  in  a clockwise  direction  (Figure  5).  The  relaxation 
factor  that  gave  the  fastest  convergence  was  1.6.  It  was  found  that  the  converged  solution 
differed  slightly  from  the  concentric  circle,  solution.  Closer  observation  of  flow  over 
symmetrical  airfoils  indicated  a dissymmetry  in  the  pressure  distribution.  For  example,  the 
location  of  the  suction  peak  was  shifted  forward  on  the  lower  surface  and  rearward  on  the 
upper  surface.  Reversing  the  direction  of  line  overrelaxation  to  counterclockwise  reversed  the 
dissymmetry.  If  the  converged  concentric  circle  solution  is  utilized  as  the  initial  condition, 
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it  satisfies  the  radial  line  version  immediately  without  further  iteration;  if  the  converged 
solution  of  a radial  line  solution  is  used  as  the  initial  condition,  it  too  satisfies  the  matrix. 
Therefore  two  significantly  different  results  would  satisfy  the  radial  line  matrix  scheme.  As 
a further  experiment,  it  was  found  that  if  a radial  line  solution  was  used  as  input  to  the 
concentric  circle  version,  several  iterations  were  required  before  convergence  was  attained. 

The  resulting  solution  was  the  same  as  that  produced  by  starting  from  incompressible 
conditions. 

Alternating  sweep  direction  from  one  iteration  to  the  next  was  also  tried.  Convergence 
was  slow  and  unsteady,  but  the  solution  was  found  to  be  more  nearly  like  that  of  the 
concentric  circle  matrix. 

This  distortion  in  the  radial  line  solution  was  considered  attributable  to  use  of  the 
updated  $ from  the  previous  line.  Use  of  the  old  value  of  $ failed  to  produce  convergence, 
but  when  the  program  was  terminated,  the  pressure  distribution  was  found  to  be  undistorted. 
As  a check  on  the  programming,  the  converged  solution  from  the  concentric  circle  matrix 
was  used  as  the  initial  condition.  Convergence  was  immediate. 

Reversed  direction  for  the  concentric  circle  method  was  also  tried.  This  consisted  of 
moving  from  the  outer  circle  to  the  center  of  the  working  plane  rather  than  the  other  way. 
More  iterations  were  required  for  convergence,  but  the  solution  was  the  same  as  for  the  other 
direction.  In  addition,  multiple  passes  through  the  matrix  solution  before  density  update 
were  tried,  but  these  offered  no  advantages  in  computational  time. 

RESULTS 

The  computer  program  of  Reference  2 was  used  to  transform  airfoils  into  the  circle. 

For  circular  and  elliptical  sections,  the  transformation  was  performed  analytically. 

VERIFICATION 

Validation  of  this  numerical  method  is  hampered  by  the  limited  number  of  accurate 
solutions  that  are  available  for  comparison  purposes.  The  critical  Mach  number  (Ml  of  a 
circular  cylinder  has  been  expressed  analytically  as  a series  summation.  Although  this  series 
has  not  been  completely  evaluated,  the  circle  was  selected  as  a check  case. 

For  the  utmost  accuracy,  a grid  of  240  x 30  was  used  and  the  free-stream  Mach  number 
was  increased  incrementally  until  local  sonic  velocity  was  reached.  Critical  Mach  number 
(0.3985)  is  in  agreement  with  the  value  reported  in  Reference  3.  The  pressure  distribution  is 
shown  in  Figure  9.  A more  practical  grid  size  of  160  x 15  gave  a result  of  0.3990. 
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Figure  9 - Pressure  Distribution  on  Circular  Cylinder 
at  the  Critical  Mach  Number  of  0.3985 

As  a check  on  the  treatment  of  circulation,  the  lift  coefficient  was  found  for  a circle  at 
which  the  two  stagnation  points  merged  and  began  to  leave  the  surface.  This  value  of  13.2 
(Figure  10)  is  in  the  range  expected  from  incompressible  flow  theory.  The  critical  value  for 
this  condition  was  used,  namely,  Mach  = 0.19.  As  a further  test,  a value  of  20.0  was  specified 
for  Cg.  As  expected,  the  stagnation  points  had  left  the  surface  and  were  found  in  the  flow 
field.  The  entire  surface  flow  was  moving  in  the  same  direction  (Figure  1 1).  No  numerical 
instabilities  were  encountered. 

Numerical  hodograph  methods  have  been  developed  that  give  a precise  (within  numerical 
error)  relationship  between  pressure  and  airfoil  geometry  for  the  design  Mach  number.  A 
design  with  a peak  Mach  number  of  0.982  was  selected  from  Reference  4.  The  agreement  is 
shown  in  Figure  1 2. 

No  lifting  hodograph  designs  were  available  for  a check  of  lifting  flow.  However,  a rough 
check  is  possible  by  determining  how  Cg  varies  with  Mach  number  for  a fixed  angle  of  attack. 
It  is  known  that  as  a first  approximation,  the  variation  is  according  to  the  Prandtl-Glauert 
factor.  The  results  for  an  NACA  0012  section  at  a = 2.0  are  shown  in  Figure  13.  The  general 
trend  matches  the  sj  1 - M2  variation. 
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Figure  10  — Mach  Number  Distribution  on  Circular  Cylinder 
for  Nearly  Critical  Flow  and  Circulation  Corresponding 
to  Merged  Stagnation  Points 


Figure  1 1 — Mach  Number  Distribution  on  Circular  Cylinder 
for  a Specified  Lift  Coefficient 
(Note  that  the  stagnation  point  has  left  the  surface) 
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LIFT  COEFFICIENT  C 


Figure  12  - Comparison  of  Finite-Difference  Solution 
with  Hodograph-Designed  Airfoil 


Figure  13  — Variation  of  Lift  Coefficient 
with  Mach  Number 
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Within  the  computer  program  itself,  lift  can  be  checked  by  comparing  the  lift  coefficient 
obtained  from  the  calculated  circulation  with  that  obtained  from  numerical  integration  of  the 
pressure  distribution.  In  nearly  all  cases,  the  difference  in  the  two  Cg  values  was  less  than 
0.2  percent  and  usually  about  0.1  percent.  This  is  an  order  of  magnitude  better  than 
obtained  by  using  the  stream  function  (Reference  1). 

COMPARISON  WITH  OTHER  METHODS 

It  is  of  considerable  interest  to  compare  results  obtained  by  compressibility  factor 
methods  with  those  of  the  present  exact  numerical  solution.  The  Prandtl-Glauert  factor  is 
given  by 

Cp  = Cp  ly/  1 - M2 

inc 

while  the  Karman-Tsien  factor  is 

Cp  = Cp  /(y/  1 - M2  +0.5  Cp  (l  -y/\  -M2  )) 

inc  inc 

where  Cp  is  the  incompressible  value  of  the  local  pressure  coefficient  and  M is  the  free- 
stream  Mach  number. 

Figure  14  compares  the  determinations  of  pressure  distribution  for  a 15-percent-thick 
ellipse  at  nearly  critical  Mach  number.  Note  that  both  compressibility  factors  overpredict 
the  nose  and  underpredict  the  midchord  suction.  Also  note  the  difference  in  the  shapes  of 
the  curves.  The  correction  factor  technique  will  always  give  a curve  which  closely  resembles 
the  shape  of  the  incompressible  distribution.  However,  the  exact  solution  will  give  a curve 
which  changes  overall  shape  rapidly  as  the  critical  Mach  number  is  approached. 

The  overprediction  of  leading  edge  suction  is  characteristic  of  the  correction  factors  and 
worsens  for  sections  with  blunt  leading  edges  that  result  in  faster  flow  acceleration.  The 
hodograph  airfoil  is  compared  with  the  Karman-Tsien  factor  in  Figure  1 5.  The  large 
discrepancy  at  the  nose  would  result  in  a very  erroneous  prediction  of  Mcr  for  this  airfoil. 
(Actually,  the  Karman-Tsien  method  would  indicate  supersonic  flow;  in  practice,  therefore, 
the  results  would  be  discarded  for  this  Mach  number.)  This  difference  is  a strong  function  of 
peak  Mach  number.  Figure  16  compares  the  Karman-Tsien  method  with  present  results  for  a 
15-percent  ellipse  with  a rounded  leading  edge  at  two  Mach  numbers.  Once  again  there  is  a 
large  discrepancy  at  high  subsonic  velocities  (Figure  16a).  For  a substantially  lower  free- 
stream  velocity,  however,  the  Karman-Tsien  approach  gives  an  acceptable  solution  (Figure  16b). 
It  should  be  pointed  out  that  the  sonic  velocity  for  lifting  cases  can  be  reached  for  free-stream 
Mach  numbers  of  0.3  or  lower. 
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Figure  14  - Comparison  of  Present  Method  with 
Prandtl-Glauert  and  Karman-Tsien  Correction  Factors 
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Figure  1 5 — Comparison  of  Karman-Tsien  Method 
with  Hodograph  Solution  of  Figure  1 2 
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Figure  16b  - At  Mach  0.400 
(Flow  condition  is  well  below  sonic  velocity) 

Figure  1 6 — Comparison  of  Karman-Tsien  Method  with  Present  Results 
for  a 15  Percent  Ellipse  with  Slightly  Rounded  Leading  Edge 


ADDITIONAL  CASES 


Solutions  by  the  present  method  for  the  NACA  0012  airfoil  are  shown  in  Figure  17  for 
different  lifting  conditions.  Since  the  0012  is  analytically  described  it  provides  a convenient 
test  case.  A circulation  control  case  for  a 1 5-percent  ellipse  and  critical  Mach  number  is 
presented  in  Figure  18.  Note  that  even  though  the  angle  of  attack  is  negative,  the  pressure 
distribution  is  such  that  positive  lift  results  from  the  specified  circulation  condition. 

CONCLUDING  REMARKS 

The  use  of  compressibility  correction  factors  for  estimating  critical  Mach  numbers  can 
lead  to  misleading  results  - even  for  comparisons  of  different  designs.  This  is  attributable  to 
the  characteristic  overprediction  at  the  leading  edge  and  underprediction  near  the  midchord 
of  suction  pressures.  Therefore,  a comparison  of  the  critical  Mach  number  of  airfoils  with 
different  types  of  pressure  distributions  could  lead  to  false  conclusions. 

The  possibility  of  extending  this  basic  method  to  supercritical  flow  was  investigated 
numerically.  The  inbedded  region  of  supersonic  flow  requires  the  use  of  one-sided  finite 
difference  equations  that  reflect  the  hyperbolic  nature  of  the  governing  equations.  (Any 
attempt  to  go  beyond  sonic  velocity  with  the  elliptic  equations  produced  severe  instabilities.) 
The  discontinuity  or  shockwave  that  terminates  the  supersonic  region  would  hopefully  be 
spread  over  two  or  more  grid  points  as  a result  of  truncation  errors,  thereby  producing  a 
continuous  flow  field  for  the  numerical  solution  to  proceed  smoothly.  The  above  concepts 
were  successfully  applied  by  Murman  and  Cole  to  the  small  disturbance  equations  (Reference  5). 

Several  changes  were  made  to  the  subsonic  program  for  this  investigation.  The  radial  line 
matrix  method  was  used  in  the  consideration  of  propagation  direction  of  transonic  disturbances. 
Rather  than  sweeping  line  by  line  in  a single  theta  direction,  the  computations  were  started  at 
the  forward  stagnation  point  and  continued  in  the  direction  of  the  flow  to  the  trailing  edge. 

At  each  point  a check  of  the  local  Mach  number  was  made;  if  it  was  greater  than  unity,  the 
hyperbolic  equations  were  used. 

No  successful  solutions  were  obtained.  Failure  occurred  shortly  after  a supersonic  region 
appeared.  The  mode  of  failure  was  that  the  Mach  number  ahead  of  the  discontinuity  grew 
steadily  to  unnatural  levels,  with  the  point  immediately  downstream  eventually  acquiring  a 
zero  or  even  a negative  velocity.  A tentative  explanation  is  that  the  truncation  error  or 
“numerical  viscosity”  in  the  equation  as  formulated  is  not  sufficient  to  smooth  out  the 
discontinuity.  However,  this  is  only  conjecture  and  no  theoretical  examination  has  been  made  of 
the  suitability  of  the  overall  concept  for  transonic  flow  solutions. 
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Figure  17b  — At  M = 0.72,  a = 0,  Cg  = 0 
Figure  17  — Present  Method  Solutions  for  NACA  0012 


APPENDIX 

FINITE  DIFFERENCE  EQUATIONS 


Finite  difference  approximations  for  the  partial  derivatives  can  be  obtained  from  Taylor 
series  expansions.  The  grid  notation  of  Figure  5 is  used  to  obtain  the  expansions  in  the 
uniform  Ar,  A0  mesh: 
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Analogous  expressions  are  obtained  for  the  r derivatives. 

The  centered  differences  used  are  given  below;  they  have  a truncation  error  of  order 
(A0)2  or  (Ar)2  as  indicated  above. 
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These  are  substituted  into  the  continuity  equation  (38).  Collecting  terms  and  rearranging 
results  in  the  linear  algeraic  equation 


where 


A$4  +B$0  +C3>2  = D 


A = p2-p4-4p0 


B = 8p0 


C = PA-P2-4p0 

i \ i 2 

D = (P4  -p2)  ( — ) [sin  (0  +a)  + Er]  + — (p{  -p3)  [r2  (<J>j  - $3) 
\ r / c 2 


- 2c  cos  ( 9 +a)]  + 2p0  rc  ($j  - 4»3)  + 4p0  r2  (4>j  + 4>3) 


The  boundary  conditions  at  the  surface  are  similarly  written  in  finite  difference  form. 
At  the  trailing  edge,  the  Kutta  condition  (0  = 0)  is 


E = 


<f>  - <T) 
‘2  *4 

2b 


- sin  (a) 


The  complexity  of  introducing  an  imaginary  externa!  grid  line  (reflection)  for  the  flow 
tangency  condition  can  be  avoided  by  using  the  following  one-sided  difference  formula: 


4>r  = cos  (0  + a)  = - — (4  4>0  - 4>3  - 3 ) 


w 


from  which  there  results  at  r - 1 


4<f>0  - $3  + 2c  cos  (0  +a) 
3 
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